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ABSTRACT 



The parameter estimation methods considered in this 
thesis are the weighted Least-Squares and Weighted Huber 
for some non-linear growth models. The properties of 
these parameter estimators derived from simulated data 
by means of (1) weighted and unweighted least-squares 
and (2) weighted and unweighted Huber robust estimation 
are compared. The error components of the simulated data 
are long-tailed and non-normal. The performance of mis-spe- 
cified models is considered. 
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I. INTRODUCTION 



The properties of the least-squares estimator for the 
linear model are well known [Ref. 2], Among its important 
properties is that it has minimum variance among the class 
of linear unbiased estimators. In Chapter II the properties 

A 

of the least-squares solution for SLnA in the growth model 

-B/t. e(t.) 

Y (t = Ae e = 1,2,. . . ,n 

will be reviewed. Then the estimate A and its statistical 
properties will be derived using the estimator for AnA by 
means of the Taylor series approximation. Thus the emphasis 
will be on estimation of the "final value" of the growth 
process . 

The least-squares estimator is easy to compute and has 
desirable properties when the assumption of normality is 
justified. However, many of the properties of the least- 
squares estimators are not robust against .non-normality as- 
sumptions (Sheffe, [3]). 

Recently (Huber, Cl]) statisticians have shown consider- 
able interest in finding estimators which are robust against 
non-normality of error assumptions. Examples of non-normal 
distributions of interest are long-tailed distributions such 
as a mixture of two normals, and also the double exponential 
and Cauchy. Several robust methods have been considered 
such as : 
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1. Using the median, 

2. Jackknifing, 

3. Trimmed mean, 

4. Huber estimator. 

Of these methods only the Huber will be considered in 
Chapter III. One major difficulty with the Huber estimation 
method is that the solution process leads to a system of 
non-linear equations which cannot, in general, be solved 
analytically. Although some asymptotic properties of the 
Huber estimator have been derived, the small sample proper- 
ties must generally be obtained by computer simulation. 

Finally, a problem is said to be misspecified if the 
data comes from the model 

- g(x, £l ) (1-1) 

but the model used is 

Y 2 = f(x,e 2 ) (1-2) 

where the functions f and g are not identical. 

Since in practice the true model is not known, it is 
important to know how good are the fits or estimates ob- 
tained by using reasonable alternative models. In Chapter 
IV the properties of several mis -specif ied model fits are 
compared when the true model is 

Bt. e(t.) 

Y( ^i) ~ A 1 + Bt. e 

l 

for A = 1 and t^ = 1,2,... 20 and B = 0.0772 using computer 
simulation. The reason for choosing B = 0.0772 is that for 
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this value of B, the largest value of E (Y) for the true model 
is 0.6, i.e. 60% of the final value (A = l) eventually reached. 
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II. LEAST-SQUARES ESTIMATION OF JlnA AND TAYLOR SERIES 
APPROXIMATION OF A 



A. LEAST-SQUARES ESTIMATION OF JlnA 
Consider the model 



Y = Ae" B/t e e 



( 2 - 1 ) 



where t is the independent variable, Y is the dependent vari- 
able, and e a random variable. 

Let the n pairs of sample observations of Y and t be 

Ctj, y x ), Ct 2 , y 2 ), ••• ,(t n , y n ). 

Assume. the hypothesis 



E(e ± ) 

E ( e i ) 



0 

O ’ 2 

e 



i - 1, 
i 7* j 
i = j 



»n 



( 2 - 2 ) 



The distribution of JlnA and the approximate distribution of 

A 

A may be found by taking logarithms of equation (2-1), i.e. 

Hn Y = Hn A - 3/t + e. * (2-3) 



The properties of JlnA can be derived by considering the 
linear model 



Z = a + 3x + u 



(2-4) 



where 

Z = Jin Y, a = Jin A, 3 = -3, x = 1/t, u = e. 

Assume the following hypothesis for the linear model 



z. = a + 3 x. + u. i = l,2,...,n 

i li > > » 
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i = 1 



j • • • y 



ECUj) = 0 
E C u i u j ) = 



0 i f j 




n 



(2-5) 



where a, 3 and cr^ are unknown parameters. The principle of 

A 

least-squares is to choose estimates a and 3 of a and 3 such 
that u| is minimized. 



The estimates a and 3 are given by 



n 



3 = 



i l 1 (x i -x) (z i -z) 



.St (X.-X ) 2 
1=1 v 1 J 



(2-6) 



where 



a = z - 3 x 



— 1 n 

x = — -2-i x. 
n i = l i 



- 1 n 
z = — • E-. z . 
n i=l l 



from these expressions it can be shown that 

A 

E (a) = a 



A 

Var (a) 




n _ 

n. E-. (x. -x) 2 

i=l v i j 



a 



2 

U ' 



(2-7) 



B. TAYLOR SERIES APPROXIMATION FOR THE LEAST-SQUARES 
ESTIMATOR A 

Let X be a random variable and Y = f(X), where it is 
assumed that the function f can be expanded in a Taylor series 
about E(X) = y; that is, 
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f(X) = f (y) + (X-y)f'(y) + 0((x-y) 2 ). 

Neglecting the higher order terms 

E(f(X)) - E( £ (y) + (X-y) £' (y) ) 

- £ (y) 

Var(£(X)) * E( (£(X) - E(f(X))) 2 ) 

- E( (£ (X) - f (y) ) 2 ) 

- E( (X-y) 2 ) • (f ' (y) ) 2 . (2-9) 

From (2-3) 

z=£nY=£nA-B/t+e 

and from (2-7) 



Var(i-nA) = 



n 

. Z 1 x? 

— lcl 1 O 2 
n e 



n. Z- (x. -x) 
i=l ^ i J 



JlnA 



Using the relationship A = e and the results of (2-8) 
and (2-9) 

y = E(X) = Jin A 
E (A) * e y 



A Q ti A A 

Var(A) = (e nA )*Var(JlnA) 



(2-10) 



n e 

n i=i f x i ' x ) 2 

Equation (2-10) shows that using the least - squares method 

A 

and the Taylor series approximation the variance of A is 
directly proportional to a 2 . 
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As an example. 



for t = 1,2,. ..,20, A = 1, and B = 10.0 
Var(A) = 0 . 0841 *Var (e) 
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III. THE HUBER ESTIMATION METHOD 



A. THE HUBER ESTIMATOR WITH SCALE UNKNOWN 

r 

Although least-squares gives the best linear (in the 
observations) estimates of the parameters in (2-3) or (2-4), 
non-linear estimates may be appropriate when the error 
terms (e*s) are long-tailed, as is often true in practice. 
Various principles for deriving estimates are possible. In 
this thesis an estimator due to Huber [Ref. 1] is utilized. 



Prob (y < Y < y + dy } 



■ 0 ( ) I dy 



C3-1) 



where Y = 0 + a e 



p(z) is the density for e. 



The likelihood function for the observations is 



n 




The log likelihood function 




n • y.; - 9 

= C log p ( — * — ) - log a ]. 



Consider 
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n 



y, - e 



J 



3L Z r P ' ^ 

3a " j=l L y, - 0 

p c -V- > 



a ) y^ - 0 

a C-l) ( - J - a — ) 3. 



To find the maximum of L where 0 and a are unknown set 3L/30 
and 3L/3a equal to zero. 



y - 0 

n P' ( 1 , ) 

j=i y7^ =0 

p ( J „ ) 



(3-2) 



in y, - 8 

- .Z. ( -3 

n j = 1 v a 



y. - 0 

P*( J a ) 

)•(-!) “ ~7i = 1 



P ( Y - 1 - 



0 



) 



In the Huber method with scale unknown, let 






(3-3) 



for some function \p to be chosen. 

The equations can be written in the form 

0 



n y, - 0 

j£i *c -4- > - 0 



, n y. 0 y. - 0 

n jfil *( -4- H -4- > - 1 



Example 1: 

Let p(z) = — i— e z ^ 

/ 2 ? 

where p is the standard normal distribution. Then 



(3-4) 



p.(z) = - z _i_ e ' zI/2 

/2tT 
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The likelihood equations can be written in the form 



n Ya 

ih - c 



0 



) = 0 



- y. - 0 y. - 0 

I n r 1 1 i r 12 

n .£ 1 1 a M a 



) = 1 



which simplifies to 

0 = — y. and 
n j=l 

a i n ^ 

a = n j-1 Cy j 0) ‘ 

These estimates are obtained by using 



iKz) = - 



* (z) = 

TTT 



(3-5) 



Thus the estimates with the normal distribution assumption 
come from the use of ty(z) = z. 

Example 2 : 

Let p be Cauchy. Then 



p(z) " ?[l l z*] and 

♦w ■ - • 

Forms of \p that are generally used 
HUBER M 

z < -ka 



iKz) = 



-ka 

z -ka < z < ka 

+ka z > ka 



(3-6) 
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SINE 



Sin( |2. ) |z| < ku 

iKz)= (3-7) 

0 otherwise. 

TUKEY 

kz(l - kz/a) | z I < ka 

I 

iKz)= (3-8) 

0 | z | > ka 

The parameter a must be selected with reference to a 
scale parameter for e. 

B. ROBUST ESTIMATION OF PARAMETERS 
Consider the linear model 

Yj = a + 8(Xj - x) + Cj j = 1 , 2 , . . . ,n 
E( £j ) - 0 

E ( e i e j ) = 0 i ? j 
E(e|) - a|. 

Let the scale factor be a. The likelihood function for the 
observations is 



L 



where x 





1 n 
- .Z, 
n j = l 



x . . 

J 




x) 



1 

• 

a 



Differentiating with respect to a, 8, and a and equating 
each derivative to zero 
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0 



3L 

3a 



n 






3L 

36 




( ^ )•(-!)•( 




0 



(3-9) 



3L 

da 




Q t y,- - a - e(x. 

C(-l)( *p )•( J 1 



a 2 



x) 




Let 

ip(z) = - - — £4 and 
yy - J p (z) 

r j = - a - 6 (x- - x) 



then the equations (3-9) can be written in the form 
n r, 



jli ’K 4 ) - 0 

'PC 4 )- (Xj - X) = 0 

jli - |] = 0. 



(3-10) 



A method of iteration for solving the system of non-linear 
equations (3-10) will now be described. First, start with 
initial values 



o(0), 6(0), a (0) . 

Re-arrange equations (3-10) : the first 

r i 

n r. n <K ) 

• ^ ( — L ) = 0 becomes .E, * r . = 0 

j=l rv a j=l r, J 

t -i’ 

and the second equation 



20 



jii *c T -i )cxj - x) = 



n PC 4 ) 

0 becomes (x. 

J = 1 r . J 



( 



) 



x ) * r j = 0 



and the third equation 



jil [iK 4- )( | J = o is now a = i ^ *( -£ )-r.. 



r i 

n *( ' aTV'i ^ — 

jh r (k) ~~ ' aCk+D-eCk+DCXj-x)] = 0 

( "4(TT } 



(3-11) 



r,(k) 

n ^( anrr ) _ — 

iT flfl— Cx j -x)[y j -aCk + l)-gCk + l)Cx.-x)] = 0 (3-12) 

C a(k) 5 



1 n r . (k+1 

a(k+l) = ± ,Z ip( -2 ) r . (k+1) . (3-13) 

n 1,-1 a (k) ^ 

From equations (3-11) and (3-12) expressions for a (k+1) and 
3 (k+1) may be obtained. These can be substituted into equa- 
tion (3-13) to obtain a(k+l). With a reasonably good guess 
for a(0), 3(0) and a(0) a few iterations should provide 
sufficiently accurate solutions of the system of equations 
(3-10) . 

Initial values for the parameters a(0), 3(0) and a(0) 
may be obtained by least squares regression. 



C. PROPERTIES OF THE HUBER ESTIMATOR 
Consider the model Y = Ae ^^e 6 ". 

Simulations were carried out with the above model for 
A = 1, B = 10, t = 1 , 2 , 3 , . . . , 20 . To investigate the properties 
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of the Huber estimation when errors are long-tailed a "wild 
shot" distribution with density function p £ (x) given by 



p e (x) - — - ■ 

e /2 tt k 



-x 2 /2k 2 , . 
e *(q) + 



-x 2 /2k 2 *9 f1 . 
e • (l _ q) 



/2tt k* 3 

was used. Four values of q were used: q = 0.9, 0.825, 0.75, 

0.6. The values of k were then chosen so that the distribution 
p £ (x) had a variance in the range 0.05 to 2.0. For each set 
of parameter values, the simulation was replicated 100 times 



to obtain 100 values of A from which the sample mean of A, 
and the sample variance were obtained. The sample distribu- 

A 

tion of A was tested for normality by means of a plotting sub- 
routine described in Appendix A. 



A 

The properties of the Huber estimate A obtained from the 
simulation are 



1. The Huber estimate A has variance which seems propor- 
tional to Var(e). Examination of Table I, Appendix B, indi- 
cates that the variance of the Huber estimator is smaller 
than the variance of the least-squares estimator. 

A 

2. The Huber estimate A has an approximate normal dis- 
tribution for 0 < Var(e)< 0.5. 

3. Fpr the "wild shot" distribution, which is symmetric, 
the Huber estimator seems unbiased. 

From Table I, Var(A) = 7.23 x 10" 3 , A = 1.0056. 

For the hypothesis E(A) = 1.0 



A-l .0 

/7 . 23x10 3 / 100 
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should have an approximate t -dis tribution with 100 degrees 
of freedom. 

The t-test at 0.10 level of significance is 



t 0 . 05 



< A - 1,0 
/7 . 2 3x10 " 3 /100 



< t 



0.95 



which may be rewritten as 



1.0 - t Q Q5 /7 . 23x10 ’ 3 / 100 < A < 1 + t Q g5 /l . 23x10 ~ 3 / 100 



t Q 95 /7.23 x10' 3 /100 = 0.0146 

0.9854 < A < 1.0146. 

“ 

The simulated value of A of 1.0056 is within the 
confidence interval identified above. 

A 

4. A good fit for the variance A against Var(c) is 
given by 

Var(A) = D Var(e) + 6*[Var(e)] 2 

where 6 is a random variable which is approximately normal 
for 0 < Var(e) < 0.5 with E(6) = 0. This is shown in Fig. 5 
where Var(A)/[Var (c) ]* is plotted against 1/Var(e). 

5.. For values of q between 0.6 and 0.9, the approximate 

' A 

relationship between Var(A) and Var(e) is 
Var(A) = 0.066 Var(e). 

This relationship is obtained from Table III, Ap- 
pendix B, by taking the arithmetic mean for q = 0.6, 0.75, 
0.825 and 0.9. The reason for using the arithmetic mean of 
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the four slopes is that they are quite close to one another 
(as compared to the spread) , and the number of replications 
(100) was not large enough to obtain a more accurate es- 
timate. In Chapter II the relationship for least-squares 
obtained from theory was 

Var(A) = 0.084 Var(e). 

Thus the Huber estimator has about 80 percent the variance 
of the least-squares model. 

6. For assumptions on the distribution of £ which 
represent a severe departure from normality, such as the 
Cauchy, examination of Table III. 2 shows that for the four 
values of the scale parameter used, the largest value of 
the variance of the Huber method for A is about 1 while that 
for the least-squares has a value of up to 10 7 . 
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IV. MISSPECIFICATION PROBLEM 



A. LEAST -SQUARES METHOD 

The theoretical values of the misspecification function 
problem when the fit is linear or log linear can be obtained 
fairly easily. First, assume that x^,X 2 , . . . ,n x are i, i,d 
with E(x) = 0 and that the function f (x^ ,X 2 , . . . > x n ) can be 
expanded in a Taylor series about the point (y^,...,y ) . 

Let fj (x 1 , . . . ,x n ) = ||- (x 1 ,...,x n ). 



Then 



n 






thus 



E[f (x 1 , . . . ,x n ) ] - £(v 1 




and 



Var(f) - E[( j S 1 (x.-y.)f.Cy 1 



9 • • • 9 






n 



- Var (x) * j 1]^ C f j (y 1 , . . . ,y n ) ] 2 . 



9 • • • 9 



Now, suppose the observed data is from the model 



Y! 

J 



1.0 t-T-4 - e j where c = 0.0772 
1 + c • j 



but the model used to fit it is 
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Then 



e. = An y. - Jin A - B/j , 
J 7 3 



and least-squares demands minimization of 

The observed values of the y^ are y! . Hence 

e.- * ln(1.0 j e e J) ■ In A ■ B/j 



- ln( ) * ej - In A * B/j. 

Defining L(A,B) as 

n 



L(A,B) = .g 1 e? 



n 



c*l 



= jSjCtoC 1+c . j ) + ej - An A + B/j] 2 . 

Differentiating this function with respect to A and B and 
setting equal to zero results in 

a n 



8L 

3L 



5 jSl [ln (T§^j) + ej - in A ♦ B/j] - 0 



55 - jll J Cln ClT^y) + e j - An A + B/j] = 0. 
The least-squares solution for An A is 



(4-2) 



AnA 



n i n r**T n i n i r • -j 

__ (ill yr ) jli^i - I7^T ] -tjli f ] ill T [ V to I§^T ] 



n n n i n i 

n t jll JT ) • c jSi y H jSi y > 

If the e are assumed to be iid with E(£j) = 0 and Var(£j) 
small, then the Taylor series approximation gives 
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E(JlnA) = 



n i n • 

( .Z.. 4 - ) r.z. J - 0 

^ 3 = 1 3 2 i ^j = l 1 +cjj J 



(4-3) 

( - Z~ i)(.Z- i Jin c* J ■ ) 

j = i iihzl j i + c»j J 



n 



n. ri + 

“( J*1 JT) - c jSj j ) 2 



For the variance of JlnA 



Mr («) ■ 



n i n l 1 

( .E- 4— ) - ( .E- i ) 

^3=1 J2- 7 v j = l]''k 



n 



1 > , 5 1 ,2 



(4-4) 



*( jfil JT)' ( jh j) 



Using (4-4) the variance for JlnA becomes 



Var(AnA)- Var(X) ^ ^ (0)) 2 

and 

Var(A) =* e^^^^ • Var (JlnA) . 



(4-5) 



1 . Comparison of Theoretical and Simulated Results 
For observed data from the model 

f 

Y j = 1,0 i+c» j c = °* 0772 

j = 1,2,. . . ,n Var(e) = 0.072 

The theoretical values are 

E (JlnA) = - 0.5768 
E(A) = 0.5617 

Var (JlnA) = 6.055 x 10 - 3 . 

The observed values for computer simulation are 

(JlnA)= - 0.580 
Var(JlnA)= 6.04 x 10 3 . 

Normality plots not included with this thesis indi- 
cated that JlnA was distributed approximately normal. Hence 
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AnA - (-0.5768) 

/O. 00604/100 

has an approximate t distribution with 100 degrees of freedom. 

The confidence interval at the 0.10 significance 
level is given by 



&nA - CO. 5768) 

0,05 ' /0. 00604/100 ' 0,95 

or 

- 0.5901 < £nA < - 0.5635. 



The observed value from simulation was £nA = - 0.580 
which is well within this confidence interval. 



B. WEIGHTED HUBER ESTIMATOR 
1. Model Y = Ae~ B/t e £ 

Quite often weights {w^ } are given to the jth obser- 
vation in order to reduce bias or the mean square error when 
it is suspected that there is misspecif ication error in the 
model. The equations (3-10) are modified to read 



jSi w j *c i- ) - 0 

n r i _ 

.Z, w. \p( -L ) (x. - x) = 0 
3=1 j ^ a J K j J 

j = l W j C ^ C )( } ' a 1 



= 0 



and the method of iteration can again be used to solve these 
systems of equations as before. 
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2. Model Y = Ae" 6 ^ e £ 

Since quite often better results can be obtained by- 

modifying the original model by the introduction of an ad- 

"B/ 1 £ 

ditional parameter the model Y = Ae ' e was modified by to 

_B / +■ c g. 

read Y = Ae e . The following is a derivation of the 

weighted Huber and Least Squares estimators for this model, 
a. Weighted Huber Derivation 

The weighted likelihood function for the obser- 
vations is 



n 

j2iC 



Any, 

C 1 



- AnA + B/t. 



) 



/• 1 -vW . 

(?) j 



For the purpose of simplifying notation replace AnA by A 
and Any^. by y^ . Then the log weighted likelihood function 
is 

n y. - A + B/t!? 

L = .£^[w. Anp( <*- ) - w. Ana]. 



Differentiating L with respect to A, B, c and a and setting 
equal to zero results in 
n 



3L 

JX 

3L 



P * 



-1 



JX ~ jil w j C “p ) C a ) = 0 



n » • i 

xx = . E. w . ( 2 — ) -±- = 0 
m 3 = 1 3 p at c 

3 -Ant 



3L _ 5 .. , p' N B 



3c 



= . E, w . ( - — ) ( 

j=l j p t^ 



1 ) = 0 



3L ,l o’ Yi - A + B / t i w n - 

3L _ t. r .. , P _ f 1 ) 1]= 0 

-a* J a 



n. t y • 

= .E-[w. ( £- ) ( -3- 

3a j=l j P 



Let = 
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Then the above equations can be written as 



n w. y. - A + B/t 



p = . E- —2- 4> ( —1 

P j=l t c 



3 



) = 0 



n w. 



q - ,-Si *r 'PC 
j x t v 



A + B/ t . 



) = 0 



(4-6) 



r = 



n w. £n t. y. - A + B/t 

- 1 r— 1 - <K 



3 = 1 



t? 

3 



i ) = o 



n 



n 



s = 



y. - A + B/t. 

a .L ± w j - j g 1 Wj (y^. - A + B/tV)i|»( -1 3_ ) = 0 



The unknown parameters A, B, c and a are the 
solutions to the system of non-linear equations (4-6) which 
can be solved by Newton's method of iteration. To simplify 
notation let r^ = y^ - A + B/t Then 

r . 



|¥- = w. ( )i|>» ( -1 ) 

9A j=l- j ^ -a J y K a J 



lE = 
3B 


n 

j-i 


w - 
J 


1 

at*: 








3 


lE = 
3c 


n 

j=i 


w. 

J 


£ ( 
a v 


12 = 


. n 


w . 


B 


3a 


iie 


J 


t C 








3 


l£ = 


n 


w . 


- ( • 


3A 


ii t 


t? 






3 




l£ = 


n 

J*1 


w . 


1 


3 B 


tj 


at 



r . 



■in t - r . 

-r— 1 )' ’ *’( ) 



■a2 



Z± 

a 



J. 

a 



r . 

-JL 

a 



r . 
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If- 



n w. R -An t. r. 

jSi ri I -c — c— ^ >»'c -i 3 

t j 3 



a<n n w i r. r -: 

3a = j = l 7^ ^' C ~i 3 

u • 

3 



3 r 

3A 



n w.Ant. n r. 

jSi -V- 1 <4) ♦*£ -i> 

t 3 



3 r 

SB 



3 r 

3c 



n w.Ant. 

jfii ^ — *■ 






at . 



c 



n w.Ant. R 

. Z 1 -J 3- - 

3 = 1 t c 

t 3 



-An t. r. 

? c — r - 1 )*'(-!•) 

£ j 



8 r 5 V nt i 

3a j = l 



t? 

3 






3_s 

3A 



n r. r. 

= -.Z- [w. (-l)iK — ■ *- ) + w.r.( — )^'( ~ ^ 

3 = 1 3 v ^ r v a ' 3 3 a 7r v a ' 



3s 

3B 



n 1 r i 1 

[w. — ip( — 1 - ) + w.r. — 

3‘1 J t? a 33 atr 
3 3 



c *’C ~i 33 



3s 

3c 



n (-Ant.) r. R (-Ant.) r. 

-. Z - [w.B 1 —'V( — •*- ) + w.r. — — - J — ip' ( —~r )] 

3 = 1 3 t c YK a J 33 ^ t c v a J 



3 



3 



is 

3a 



jil w j - jil ( r^) *’C r ~i 



T T 

Let z = (z 1 > z 2 ,z 5 > z /[ ) = (p,q,r,s) 



T T 

x = (x 1 ,x 2 ,x 3 ,x 4 ) = (A,B,c,a). . 
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Then the iterations are given by 



Vi - A i * AA i 



B i+1 ' B i * AB i 





C i + 1 


= c 


. + 
1 




a i + l 


= a 


♦ 

1 


where 


the 


A ' s 


are 


9z j 


9z j 


9z^ 


9z 


9x^ 


9x 2 


8X 3 


3x 


9z 2 


9z 2 


9z 2 


3 z 


9x^ 


9x 2 


8x 3 


3x 


9z 3 


9Zj 


9Zj 


3z 


9x^ 


9x 2 


8X 3 


3x 


9z 4 


8z 4 


8z 4 


3z 


9x^ 


9x 2 


9x 3 


3x 













AA. 

l 




p ! 




AB. 

l 




q 




Ac^ 




r 




Aa . 




s 




l 












_ 



b. Weighted Least-Squares Derivation 

The weighted least-squares estimator for the model 

Y = Ae" B//t e £ 

lead to' the normal equations 

n g 

p = . E- w.[£ny. - AnA + — — ] 

F 3=1 3 J t c 

3 



= 0 



q = w. -p [Any. - UnA + -p ] = 0 
J ± J tV J tT 
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r = . w. — £■ £nt.[£ny. - £nA + — |r ] = 0 
^ ^ t • ^ t • 



J 



and the iterations are given by 

£ny : 



c . . = c . 
l + l 1 



[ n 
L £ 

_ j = l 



J. 



Vi 



(JlnA) 



n 



i . t . c . 
j=l j l 



w i 

1 — + B . 5 



w . 



i . t . 2c . 
j=l j l 



C - B. 



Ant . 

W. T — L ] 



i j=l j tj Ci 



n w.£ny. 
C.S- l 



Vi 



B. - 
1 



j* 1 Vi*! 



(In A) 



w . 

J. 



n 



j' 1 



+ B . . E 



w . 

JL 



i j=l t.2c. 



i+1 



(XnAj i 



n 

.2 



w . 



j =1 (Vi+l 5 



(AnA) i + 1 =CJlnA) i - 



n w. Any. Ant. 
[ .Z- -J 3 1 

l J=l tjC. 



- (AnA) 



w . ilnt . 

J I 



'j i + 1 



1 J=1 Vi + i 



B 



n 

.Z. 



w 



i+1 j=l tj2 



n 



w . ilnt . 

[ j-^l 



'j i + 1 



Bt 



C. COMPUTER SIMULATION OF THE MISSPECIFIED MODELS 

The data analyzed was generated from the model Y = 1 *° 1+Bt € 
where B = 0.0772 with var(e) = 0.1, 0.2 and 0.3 initially. The 
value of B was chosen so that the expected values have a range 
from 0.07 to 0.60. The values for t were from one to twenty. 

"a a 

The sample mean A and sample variance Var(A) were calculated 
for a sample size of 100. The weighting system chosen was 
w(j) = 1.0, j, j 1 * 5 , j 2,0 , 0. , O.SS^ 20- ^ in order to 
find weights which show promising results so as to concentrate 
the analysis on the interesting cases. 
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1. Numerical Results 



that 



Preliminary examination of Tables VI, IX, and X show 



- B / t^ e 

1. For the least-squares fit to Y = Ae e the 



expected value of A is practically independent of the system 

A 

of weights chosen with E(A) approximately equal to 1.30, 
which is quite surprising. 

“ C / G 

2. For the same model Y = Ae ' e the least-.squares 
method and the Huber method have approximately the same ex- 
pected value and variance for low weights. 

3. For heavy weighting systems, the bias of the 
Huber method is considerably less than the bias of the least- 
squares method. 

4. The iteration for the Huber method with the model 

_B/t C £ 

Y = Ae ' e for varCc) = 0*2 had some peculiar features. 
Specifically, 



a. About eighty to ninety percent of the simu- 
lations with different random number seeds did not result in 
convergence . 

b. Examination of Table V appears to show that 
the system of non-linear equations for the Huber method have 
multiple solutions. 

c. When the initial starting values for the 
Newton iteration are the convergent values of one of the 
sequences, but if a different seed is used, divergence of 
the iteration still usually occurred. 
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The probable cause of the divergence for the itera- 
tion of the Huber method for the mis - specif ied model is that 
the var(e) chosen was too large. One of the difficulties of 
the Huber method for the model Y = Ae ' e is chosing a 
suitable set of parameters (4) as the initial starting values 
for the iteration. Usually it is quite difficult to chose 
even three suitable initial values. 

As a result, it was decided to concentrate further 
analysis of the mis -specif ied problem using heaving weighting 
systems with the Huber method for the model Y = Ae ' e as 
these gave the least bias. Also, as the initial value of 
var(e) used was too large to be practical, the values of 
var(e) was reduced to var(e) = 0.005, 0.01, and 0.02. The 
weighting systems chosen were w(jj = j , j , j , j , o.85 v JJ , 
0 . 7 < 20 -«. 

Examination of Tables VII and VIII shows that the sys- 
tem which weighted the later data values heavily and early 
data values very lightly still gave the best results. The 
weighting system given by w(j) = has the smallest bias and 

A 

mean square error although it had an estimator A with the 
largest variance. 

2 . Individual Curve Fitting Trials 

The following curve fitting experiments were performed 
for the Huber method and the least-squares method for the 

“ C / 1 £ 

model Y = Ae e with data generated from the model Y = 1.0 
(B*t/1+B*t)e with parameter values previously mentioned. Two 
curves were plotted in each figure. The curve which has the 
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smaller value of y for small values of t is the fitted curve. 
The second curve is the expected value of T for the data 
generated. The figures are: 

1. 6-13 var(e) = 0.005 

2. 14-21 var(e) = 0.010. 

The unweighted least-squares curves (Figures 6,14) 
were obviously bad fits. Comparison of the various figures 
show that in general the heavier weighting systems gave 
better fitted curves for large values of t and poorer fits 
for small values of t. 
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V. DISCUSSION 



A. JUSTIFICATION FOR HUBER METHOD 

For the currently specified model given by Y = Ae ^^e 6 " 
with E(e) =0, the Huber method has a slightly smaller vari- 
ance than least-squares for the random variable e consisting 
of a mixture of two normals. For Var(e) < 0.5, the variance 

A 

of the Huber estimator A is approximately normal. Thus the 
Huber estimator of A is a better estimator than least-squares 
for long-tailed distributions. For non-symmetric distribu- 
tions of e the least-squares and Huber methods are compar- 
able (see Table IV). 

Another possible objection to the least-squares method 
for the fit Y = Ae ' e is that very often it is not known 
what the true model is. The least- squares method does not 
seem robust against a mis -specified model, judging from our 
experiments. Thus the Huber method might well be preferred 
to the least-squares method when the true model is not known, 
judging from the experiments carried out to date. 

— ^ J ^ 0 

Another variation of the model Y = Ae ' e is Y = A* 
e ^^ t+ ^e e . For the data Y = 1.0 (Bt/1+Bt)e e was considered. 
This model did not appear promising for the following rea- 
sons : 

1. Comparing the plots for the various methods of fit 
* B / 1 0 

for Y = Ae ' e shows that the fitted curves were too flat 
as compared to the curve Y = 1.0 (Bt/l+Bt) for values of t 
satisfying 12 < t < 20. Since the parameter k was introduced 
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to match more closely the slopes of the fitted curve and the 
data for values of t between say 15 and 20 as well as the y 
values, this implies that k would have to be negative. How- 
ever, the function y = Ae has an infinite discontinuity 

at t- = -k. Thus values of k < -1 are not reasonable for the 
range of values of t considered above. One way of overcoming 
this difficulty would be to discard the observations for 
0 < t < 5 , since the future values are less dependent on the 
distant past observations than those observations which are 
more recent. 

B. AREAS OF FURTHER STUDY 

t 

1. Since the exact form of the growth is not known, 

further investigation in this area. For example, the model 
used to generate the data in this thesis was Y = A* e • 

— Ct £ 

However, there is no reason to believe that Y = A(l-e )e 
or any similar growth model could not have been used. 

2. The range of values for t was from one to twenty 
(1-20). Other ranges on t need to be investigated. 
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APPENDIX A: 



SUBROUTINE NMPLOT 



The subroutine NMPLOT was used as an indicator as to 

whether a set of n observations x, ,x 0 , . . . ,'x come from a 

12 n 

normal distribution. To carry out the test - , the numbers 

x^,X 2 »...,x n are sorted in increasing order say y^ ^ , . . . ,y . 

— 1 n 

The sample mean y = — .£,y. anc ^ t * ie sam Pl e variance _ 

, i n _ y i - y 

s - - — r (y- - y) are then calculated. Let z- = — * 

n-1 j =1 w j ' ' j s 

If the original set x^,...,x n were normal, then for rather 
large n (at least 50) the set z . , j=l,...,n should have 
properties similar to that of the ordered statistics from a 
standard normal distribution. Let $ be the cumulative dis- 
tribution function from the standard normal. Then a plot 
of j/n versus $(Zj) should lie approximately on a straight 
line passing through the origin with slope one. To decide 
whether the sample x^,...,x is approximately normal, a plot 
obtained from n observations from a standard normal and the 
two plots compared. The subroutine is part of the computer 
program for the weighted Huber method for the model Y = Ae e . 



39 



Wild Shot Distribution 
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Wild Shot Distribution q = 0.75 
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Table II. Comparison of Least-Squares with Huber Method q = 0.75. 



Var (e) 


q = 0 . 6 
Var ( A ) 


0.1 


0.0060 


0.2 


0.0144 


0.3 


0.0257 


0.4 


0.0222 


0.5 


0.0327 


0.6 


0.0393 


0.7 


0.0458 


0.8 


0.0492 


0.9 


0.0668 


1.0 


0.0512 


1.1 


0.0697 


1.2 


0.0871 


1.3 


0.0807 


1.4 


0.0752 


1.5 


0.1381 


1.6 


0.1091 


1.7 


0.0933 


1.8 


0.0982 


1.9 


0.1493 


2.0 


0.1685 


Slope 


0.0641 



Y = 1.0 e' 10/t e e 



q=0. 75 
Var (A) 


q=0 .825 
Var (A) 


0.0058 


0.0072 


0.0118 


0.0145 


0.0254 


0.0235 ' 


0.0197 


0.0250 


0.0181 


0.0255 


0.0273 


0.0345 


0.0489 


0.0566 


0.0468 


0.0451 


0.0515 


0.0679 


0.0876 


0.0734 


0.0552 


0.0802 


0.0427 


0.0614 


0.0588 


0.0911 


0.1080 


0.1017 


0.0789 


0.0888 


0.1060 


0.0858 


0.0860 


0.1008 


0.1700 


0.1071 


0.1330 


0.1270 


0.1730 


0.2015 


0.0593 


0.071 



Table III.l. Variance for the Huber Method 
Distribution . 



q=0 . 9 
Var (A) 

0.0068 

0.0156 

0.0209 

0.0241 

0.0329 

0.0434 

0.0490 

0.0508 

0.0553 

0.0683 

0.0681 

0.0939 

0.0909 

0.1340 

0.1530 

0.1389 

0.1214 

0.1970 

0.2070 

0.1650 

0.070 

Wild Shot 
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e is mixture of two distributions 

Z with probability 0.9 

Z^ with probability 0.1 

X is the standard normal 

HUBER 



LEAST SQUARES 



Z 1 


P 

0.0707 


A 

1.020 


Var(A) 

0.00100 


A 

1.021 


Var (A) 
0.00099 


0.25 + pX 


0.1 


1.029 


0.00147 


1.030 


0.00147 




0.1414 


1.030 


0.00218 


1.032 


0.00215 




0.0707 


1.023 


0.00102 


1.024 


0.00101 


Z 2 =0.25 


0.10 


1.030 


0.00142 


1.030 


0.00141 




0.1414 


1.020 


0.00203 


1.020 


0.00201 



Table IV. Comparison of Least Squares and Huber Methods for 
Non -Symmetric Distributions. 
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Table V. 


Convergent Sequences for the Huber 
Mis-specif ied Model. 

Y = Ae e e with weights w(i) = 


Method 

1 


Iteration 


(£nA) ± 


B. 

l 


C i 


a. 

i 


1 

• 


0.2000 


3.000 


0.5000 


0.8000 


• 

16 


-0.1083 


3.6682 


0.7245 


-18.179 


17 


-0.1086 


3.6678 


0.7246 


-18.072 


18 


-0.1086 


3.6678 


0.7246 


-17.965 


19 


-0.1086 


3.6678 


0.7246 


-17.858 


15 


0.7872 


3.1258 


0.2924 


9.006 


16 


0.7969 


3.1349 


0.2911 


8.226 


17 


0.7971 


3.1351 


0.2910 


7.437 


18 


0.7971 


3.1351 


0.2911 


6.649 


19 


0.7971 


3.1351 


0.2911 


5.862 


20 

* 


0.7971 


3.1351 


0.2911 


5.076 


* 

15 


0.0912 


2.6268 


0.4214 


0.8348 


16 


0.0912 


2.6268 


0.4214 


0.7991 


17 


0.0912 


2.6268 


0.4214 


0.7657 


18 


0.0912 


2.6268 


0.4214 


0.7348 


19 . 


0.0912 


2.6268 


0.4214 


0.7064 


20 

• 


0.0912 


2.6268 


0.4212 


0.6803 


• 

15 


0.3407 


3.2417 


0.4544 


0. 1975 


16 


0.3346 


3.2353 


0.4560 


0.1929 


17 


0.3360 


3.2367 ■ 


0.4557 


0.1939 


18 


0.3356 


3.2363 


0.4558 


0.1936 


19 


0.3357 


3.2364 


0.4557 


0.1937 


20 


0.3357 


3.2364 

45 


0.4558 


0. 1936 





Y = 1.0 


Bt 


e £ b - 


0.0772 q = 


0.9 




1 + Bt 




w(j) 


Var (e) 




Var (A) 


T 

A 


M.S.E. 


1 


0.10 




0.00041 


0.569 


0.186 


3 


0.10 




0.0055 


0.681 


0.107 


j 1 ' 5 


0.10 




0.0123 


0.697 


0.104 


j 2 


0.10 




0.0220 


0.740 


0.090 


0 . 70 20 


0.10 




0.0414 


0.761 


0.098 


0.85 2 °-5 


0.10 




0.0038 


0.655 


0.123 


1 


0.20 




0.0154 


0.560 


0.195 


3 


0.20 




0.0149 


0.663 


0.128 


j 1 - 5 


0.20 




0.0325 


0.722 


0.110 


3 2 


0.20 




0.068 


0.78 


0.116 


Q.70 20 ~i 


0.20 




0.127 


0.839 


0.153 


o.ss 20 '^ 


0.20 




0.0138 


0.674 


0.120 


1 


0.30 




0.0209 


0.570 


0.206 


3 


0.30 




0.0183 


0.645 


0.144 


j 1 . 5 


0.30 




0.0344 


0.701 


0.124 


.2 

1 


0.30 




0.066 


0.754 


0.127 


o. yo 20 "^ 


0.30 




0.199 


0.852 


0.221 


o.ss 20 "^ 


0.30 




0.0134 


0.663 


0.127 



Table VI. Weighted Huber Method with Mis -specified Model 
Y = Ae‘ c/t e t . 
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q = 0.9 



Y = 1.0 e £ B = 0.0772 



w(j) 


Var (e) 


Var (A) 


A 


M.S.E. 


O^O 20 ' 3 


0.005 


0.00333 


0.757 


0.0629 


o.ss 20-3 


0.005 


0.00043 


0.662 


0.115 


.1.5 


0.005 


0.00056 


0.710 


0.0847 


i 2 


0.005 


0.00118 


0.749 


0.0642 


o 

• 

-*4 

O 

ts> 

O 

l 

«— i. 


0.01 


0.00503 


0.756 


0.0646 


0 . 85 20 


0.01 


0.00087 


0.664 


0.114 


j 1 . 5 


0.01 


0.00121 


0.718 


0.0807 


j 2 


0.01 


0.00154 


0.750 


0.0640 


0. 70 20- ^ 


0.02 


0.0119 


0.781 


0.0600 


0 . 8 S 2 0 ~ 


0.02 


0.00117 


0.659 


0.117 


.1.5 


0.02 


0.00253 


0.707 


0.0883 


j 2 


0.02 


0.00432 


0.742 


0.0709 



Table VII. Weighted Huber Method with Mis-specified Model 
Y = Ae‘ c/t e e . 
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Y - 1.0 e £ B = 0.0772 q = 0.9 



a A a 



w(j) 


Var (e) 


Var (A) 


A 


M.S.E. 


j 3 


0.005 


0.0033 


0.795 


0.0453 


j 4 


0.005 


0.0043 


0.835 


0.0315 


j j 3 


0.01 


0.0035 


0.796 


0.0451 


; 4 


0.01 


0.0077 


0.816 


0.0416 


j 3 


0.02 


0.0092 


0.790 


0.0533 


j 4 


0.02 


0.0179 


0.839 


0.0438 



Table VIII. Weighted Huber Method with Mis-specif ied Model 
Y = Ae' c/t e £ 
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Y = 1.0 


Bt 


e £ B = 


0.0772 q = 


0.9 




1 + Bt 


w(j) 


Var (e) 




Var (A) 


A 

A 


M.S.E. 


1 


0.10 




0.00430 


1.30 


0.094 


j 


0.10 




0.00415 


1. 33 


0.113 


J 1 - 5 


0.10 




0.0067 


1.30 


0.097 


.2 

J 


0. 10 




0.0074 


1.29 


0.092 


0.70 Z0 ~j 


0.10 




0.00683 


1.26 


0.074 


0.85 20 'j 


0.10 




0.0030 


1.30 


0.093 


1 


0.20 




0.00896 


1.28 


0.087 


j 


0.20 




0.0126 


1.34 


0.128 


j 1 . 5 


0.20 




0.0146 


1.30 


0.105 


.2 

J 


0.20 




0.0182 


1. 32 


0.121 


0.7 20 -j 


0.20 




0.0154 


1.25 


0.078 


O.SS 20 ^ 


0.20 




0.0140 


1.32 


0.116 


1 


0.30 




0.0165 


1.28 


0.095 


j 


0.30 




0.0130 


1.34 


0.129 


j 1 . 5 


0.30 




0.0179 


1.31 


0.114 


.2 


0.30 




0.0268 


1.29 


0.111 


0. 70 20_;) 


0.30 




0.0361 


1.27 


0.109 


0 . 85 20 


0.30 




0.0139 


1.31 


0.110 



Table IX. Least -Squares Method with Mis^pecified Model 
Y = Ae _D / tC e e 
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B = 0.0772 



q 



0.9 



Y = 



1.0 



Bt 

1+Bt 



A A A 



w(j) 


Var (e) 


Var (A) 


A 


M.S.E. 


1 


0. 10 


0.00285 


0.566 


0.191 


3 


0.10 


0.00526 


0 . 654 * 


0.125 


.1.5 


0.10 


0.00685 


0.673 


0.114 


3 2 


0.10 


0.0152 


0.738 


0.084 


0.70 20 ~i 


0.10 


0.0194 


0.755 


0.079 


0 . 85 20 -J 


0.10 


0.00677 


0.652 


0.129 


1 


0.20 


0.00517 


0.572 


0.188 


3 


0.20 


0.0116 


0.654 


0.131 


j 1 . 5 


0.20 


0.0214 


0.711 


0.105 


.2 

3 


0.20 


0.0200 


0.694 


0.114 


•r-» 

1 

O 

cvi 

r^. 

• 

o 


0.20 


0.0599 


0.741 


0.127 


0 . 85 20 'j 


0.20 


0.0139 


0.658 


0.131 


1 


0.30 


0.00831 


0.570 


0.193 


3 


0.30 


0.0165 


0.668 


0.127 


. 1.5 


0.30 


0.0181 


0.673 


0.125 


.2 

3 


0.30 


0.0429 


0.751 


0.105 


O 

• 

o 

ts) 

0 

1 

u. 


0.30 


0.0597 


0.733 


0.131 


0 . 85 20 'j 


0.30 


0.0230 


0.658 


0.134 


Table X. 


Least-Squares 


Method with 


Mis-speci£ied 


Model 



Y = Ae' c/t e e 
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Figure 1. Comparison of Variances of Huber and Least-Squares 
Estimator, var(e) < 0.5. 

V 



i 
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Figure 2. Linear Fit for Variance of Huber Estimator for 
q = 0.9. 
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Figure 3. Linear Fit for Variance of Huber Estimator for 
q = 0.825. 
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Figure 4. Linear Fit for Variance of Huber Estimator, for 
q = 0 . 60 . 
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Figure S. Linear Fit for Variance of Huber Estimator with Transformed Variables 

n = H Q * 
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I 



Figure 6. Unweighted Least-Squares, Var(e) = 0.005. 
A = 0.572. 
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Figure 7. Unweighted Huber, Var(e) = 0.005, 
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>> 
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Figure 10. Weighted Least-Squares, var(e) = 0.005, w(j)»j 3 
A = 0.797. 
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/ 



Figure 11. Weighted Huber, var(e) = 0.005, w(j) = j 3 
A = 0.810. 
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Figure 12. Weighted Least-Squares, var(e) = 0.005 , 

A = 0.824. 

* * 

\ 
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Figure 14 . Unweighted Least-Squares, var(e) = 0.10, A = 0.548. 



i 
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Figure 15 .Unweighted Huber, var(e) = 0.010, A = 0.548. 
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Figure 16. Weighted Leas t- Squares , var(e) = 0.010, w(j) 
= j , A = 0.664. 
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Figure 18. Weighted Least-Squares, var(e) = 0.010, w(j) = 
j 2 , A = 0.695. 







68 




69 




;Figure 20. Weighted Least-Squares, var(e) = 0.010, w(j) 
= y , A * 0.691. 
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Figure 21. Weighted Huber, var(e) = 0.010, w(j) = j 4 
A = 0.691. 
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( ( I ) 10 HSM)dX 3 *( ( ( I ) 1 * 30 * 10 + 0 * 1 )/( ( 1 ) 1 * 30 * 10 ) ) 



Y ( I ) = ALOG( YY ( I ) ) 
xin = -(l.o/Tun 
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WEIGHTED HUBER METHOD FOR THE MODEL Y=A*EXP ( -B/T ) *EXP ( U ) 
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